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Abstract 

We present a non-iterative method to deconvolve the spatial response function or 
the point spread function (PSF) from images taken with the Suzaku X-ray Imaging 
Spectrometer (XIS). The method is optimized for analyses of extended sources with 
high photon statistics. Suzaku has four XIS detectors each with its own X-ray CCD 
and X-Ray Telescope (XRT) and has been providing unique opportunities in spatially- 
resolved spectroscopic analyses of extended objects. The detectors, however, suffer 
from broad and position-dependent PSFs with their typical half-power density (HPD) 
of about 110". In the authors' view, this shortcoming has been preventing the high 
collecting area and high spectral resolution of Suzaku to be fully exploited. The 
present method is intended to recover spatial resolution to ~ 15" over a dynamic 
range around 1:100 in the brightness without assuming any source model. Our de- 
convolution proceeds in two steps: An XIS image is multiplied with the inverse re- 
sponse matrix calculated from its PSF after rebinning CCD pixels to larger-size tiles 
(typically 6" x 6"); The inverted image is then adaptively smoothed to obtain the 
final deconvolved image. The PSF is modeled on a ray-tracing program and an ob- 
served point-source image. The deconvolution method has been applied to images of 
Centaurus A, PSR B1509-58 and RCW 89 taken by one XIS (XIS-1). The results 
have been compared with images obtained with Chandra to conclude that the spatial 
resolution has been recovered to ~ 20" down to regions where surface brightness is 
about 1:50 of the brightest tile in the image. We believe the spatial resolution and 
the dynamic range can be improved in the future with higher fidelity PSF modeling 
and higher precision pointing information. 

1 



Key words: method: data analysis - - technique: image processing X-ray: 
general 

1. Introduction 

Three X-ray astronomical observatories, Chandra, XMM-Newton, and Suzaku, are oper- 
ational in orbit now. All have one or more X-ray imaging spectrometers, each made of an X-ray 
mirror and an X-ray CCD array Among the three, Suzaku has the second largest effective area 
at higher energies (590 cm 2 at 8 keV) when its four X-ray imaging spectrometers (XISs) are 
combined, the best energy resolution (FWHM ~ 130 eV at 5.9 keV) and the lowest background 
rate 1.0 x 10 -7 counts s _1 keV -1 arcmin -2 cm -2 at 6 keV) (Mitsuda et al. 2007). However, 
these advantages are often compromised by its relatively poor spatial resolution. The half power 
diameter (HPD) of the Chandra, XMM-Newton and Suzaku X-ray imaging spectrometers are 
about 1", 8", and 110 /y respectively. The point spread function (PSF) of the four Suzaku XISs 
is not only large and complex but depends on the position in the focal plane making analyses of 
source-rich regions or extended sources difficult (Serlemitsos et al. 2007). Whereas the HPD 
of the Suzaku X-ray telescopes (XRTs) is large, their PSF has a sharp core with an exponential 
radial profile of characteristic spread ~ 10" (Serlemitsos et al. 2007). The pixel size of the XIS 
CCDs is l // .04 and much smaller than the width of the PSF (Koyama et al. 2007a). If a high 
number of photons are available, spatial resolution of XIS can be recovered to the PSF core 
size {p=L 10") . We present here a first attempt to improve spatial resolution of Suzaku XIS by 
deconvolving the PSF. 

In astronomy, many attempts have been made in the past to obtain better spatial in- 
formation from observed data as reviewed by Starck et al. (2002) and Puetter et al. (2005). 
Several attempts have been reported for ASCA X-ray images which had a PSF similar to that of 
Suzaku. White & Buote (2000) employed a maximum-likelihood method to deconvolve images 
of the Gas Imaging Spectrometer (GIS) aboard ASCA. They reconstructed spatially resolved 
spectrum by first deconvolving energy-selected images with a maximum-likelihood method and 
then reassigning individual observed photons to a position in the deconvolved image space 
with a Monte-Carlo method. Hwang & Gotthelf (1997) applied a Richardson-Lucy method 
to Solid-state Imaging Spectrometer (SIS) images of a supernova remnant in selected energy 
bands. In these methods, images are deconvolved by assuming a model of extended emission 
or a collection of point sources. 

Because of the high through-put provided by the four Suzaku XISs, photon statistics 
does not limit fidelity of image deconvolution for bright targets. The inverse matrix method 
used in this work can recover the true image if the response matrix is known accurately and if 
the photon statistics is high in the region of interest. The method has a merit that it does not 
require any prior model to fit and hence relatively free of systematic bias. It is best suited for 
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complex extended sources. 

Our prime targets are galaxy clusters, pulsar wind nebulae and supernova remnants. 
The response function of the detector must be modeled accurately and the noise due to Poisson 
fluctuation must be controlled well to reproduce the image faithfully in typical brightness 
variation (an order of magnitude) in these extended objects. Since the PSF of Suzaku XIS 
is complex, the method requires intense labor. Once a procedure is established, however, this 
method can be set up for automatic deconvolution with minimum human intervention. 

The goal of this deconvolution is to recover spatial resolution to ~ 15" while keeping 
fidelity in a dynamic range around 1:100 in the brightness. Deconvolution proceeds in two 
steps: Each XIS image is multiplied with the inverse response matrix calculated from its PSF 
after rebinning the raw CCD image to larger tiles (typically 6" x 6"). This is required to secure 
high photon counts in each image elements. The inverted image is then adaptively smoothed 
to obtain the final deconvolved image. 

When this work began, the released XIS images were degraded due to wobbling of 
satellite pointing (Serlemitsos et al. 2007). We have added one more step to correct for this 
pointing error by using a bright point source. We have modeled the PSF in two ways: one 
based on the ray-tracing simulation program developed by the Suzaku team (Ishisaki et al. 
2007) and the other by fitting the image of Centaurus A. The two PSFs agree well for XIS-1 
but not for the others. For this reason we use only one of the 4 XISs when deconvolving the 
images of Centaurus A, PSR B1509-58 and RCW 89. 

The paper proceeds as follows. We present our image-deconvolution method in section 2. 
The Suzaku data processing and the correction prodedure for the pointing error are described 
in section 3. The PSF modeling is explained in section 4. We then apply our deconvolution 
method to XIS-1 images and compare the results with corresponding images of Chandra ACIS 
in section 5. Conclusions and our future plan are given in section 6. 

2. Deconvolution by Inverse Response Matrix and Adaptive Smoothing 

We here denote an image derived by multiplying an observed image with the inverse 
response matrix as an "inverted image". When an "inverted image" is adaptively smoothed, 
we called it an "deconvolved image." 

2.1. Inverse Matrix Method 

The relation between an observed image d(x) and a true image on the sky s(x) is 
represented by 



where P(r'\r) is PSF at the source position r' and n{f) represents collection of "noises" like 
Poisson fluctuation in photon counts, instrumental noise, and imperfect PSF modeling. 




(i) 
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An observed image taken by a pixelized imager such as a CCD is represented by a finite 
dimension vector d = {di}i=i..jv, where TV is the number of pixels used in the inversion. We 
combined multiple CCD pixels into a tile to control Poisson fluctuation. The image region 
contains 64 x 64 tiles and hence TV = 4096. 

The discretized version of equation (1) is 

di = ^PijSj + Ui (i = 1...JV), (2) 
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where pj = {pij}i=i...N is the PSF at the j-th tile, fj. The response matrix, is then a matrix 
consisting of N PSF- vectors and equation (2) becomes 

d = Ps + n, (3) 

where P is, 

P = {PiP2-Pn} (4) 

The inverse response matrix is represented by P _1 . The inverted image s' = {sj}i=i...jv is 
calculated by multiplying the inverse response matrix P _1 (dimension TV x N) to a raw image 
vector d as, 

s' = P~ l d = P-\Ps + n) = s + P~ l n (5) 

When the noise term P _1 n is negligible relative to the signal s, the inverted image is 
a good approximation to the true image. However, this is not true generally. In many cases, 
the noise term, P _1 n, dominates over the signal in the inverted image as described in Starck 
et al. 2002; Puetter et al. 2005. We control this noise term by adaptive smoothing as will be 
described in the next subsection. 

2.2. Adaptive Smoothing of Inverted Image 

We employ a technique known as adaptive smoothing to control the noise term in 
the inverted image while keeping highest spatial resolution compatible with photon statistics. 
Adaptive smoothing is a generic smoothing method where the spatial resolution is balanced to 
the signal-to-noise ratio expected for each position in the image (e.g. Lorenz et al. 1993; Huang 
& Sarazin 1996; Ebeling et al. 2006). It can be optimized to a strategy: we take a uniform 
significance approach in which the smoothing scale is adjusted so that the smoothed data have 
a similar signal-to-noise ratio everywhere in the image. This approach has been employed in 
AKIS (Huang & Sarazin 1996) and ASMOOTH (Ebeling et al. 2006). In the present adaptation, 
we smooth the inverted image not the observed image. 

Our adaptive smoothing is performed by multiplying a multi-scale smoothing kernel, 
/C(a, r) tile-by-tile over the inverted image plane. The relation between the input image s'(r) 
and the smoothed image s*(r) is represented by 

s*(r) = y 4 5 / (r / )/C(a(r);r / — r)dr' (6) 
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We here note that the width a(r) is defined in the smoothed image space r, not in the input 
image space r' . 

We use a Gaussian kernel. The smoothing matrix (TV x TV dimension) for pixelized image 

is 

F = {/C(a(r;),r - - u)} l=1 Nj=1 N (7) 

i=l...NJ=l...N 

The smoothed image s* and the noise image n* are represented by 

£ = FP~ 1 d = Td (9) 

n* = FP- 1 n = Tn (10) 

where FP~ l — T— {ty}i=i...jvj=i...JV- The signal-to-noise ratio at each pixel {STVit^j^. jv is 
then, 

5 W = 4 = ^=§iM (li) 

When Poisson noise is dominant and <^ is large, = y/dl. Since the Poisson noise of photon 
counts at each pixel can be regarded as mutually independent, the equation (11) is reduced to 
the following equation. 

SNRi = E ^ 2 = ^Ml (12) 
Search for Optimal Smoothing Scales 

Our strategy is to find a set of smoothing widths {o"i}i=i..jv so that the signal-to-noise 
ratio {SNRi} i= i„N agree with a given preset value SNR opt at every pixel on the image. We 
developed a deconvolution program to search for an optimum iteratively. A similar method has 
been used in ASMOOTH (Ebeling et al. 2006). The program starts to test smoothing widths 
{ci}i=i...N from an initial minimum value. The larger the smoothing width is, the lower the 
spatial resolution of the inverted image will be. The program increases each of {(Ji}i=i...N step 
by step, then stops to increase Oi when SNRi reaches the preset desired value. This step is 
repeated until SNRi reaches the desired value for all pixels or Oi reaches a preset maximum 
value. 

The deconvolution program takes four parameters: the desired signal-to-noise ratio 
SNR opt , the minimum smoothing scale a m i n , the maximum smoothing scale a max , and the 
increment for each iterative step Act. Among the four parameters, the deconvolution result is 
sensitive only to the signal-to-noise ratio SNR opt . The minimum scale a m i n and the increment 
for each step Aa are chosen to be smaller than 1 tile size of the input image (= 6") to achieve 
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the maximum resolution. The maximum scale a max is chosen as large as the standard devia- 
tion of the PSF (~ 30"). We use the parameter values, SNR opt = 4.0, a min = 2".4, a max = 24", 
Act = 1".4 in the tests of the deconvolution method in section 5. 

3. XIS Data 

We used Suzaku archival data released via the pipeline processing version 1.2.2.3. 
The data has been processed using the HE AD AS software version 6.1.2 released from the 
NASA/GSFC Guest Observer Facility. We extracted raw XIS images from the XIS event data 
screened by the standard screening procedure. 

3.1. Correction for XRT Alignment Error 

The pointing direction of the Suzaku's XRTs has been found to wobble by ~ 40" roughly 
synchronized with the 96-minute orbital motion of the spacecraft (Serlemitsos et al. 2007). This 
wobbling is now understood as due to thermal distortion. The distortion is introduced when 
the side panel #7 on which the start trackers and the gyroscopes are mounted, is illuminated 
by sun-lit Earth. (Serlemitsos et al. 2007). Software to correct for the alignment error is now 
available (Uchiyama et al. 2008). However, it was not available at the time of this work and 
we corrected for this error using X-ray point sources in the data. 

The corrections are done in two steps. Firstly, data was divided into 500-second intervals 
and the position of a bright point source is monitored for each interval. Figure 1 shows the 
time variation of the source position in the observed Cen A data in the initial performance 
verification phase. The periodicity of the 96-minute period is seen clearly. The images for 
500-second intervals were shifted so that the Cen A align to one position. Figure 2 shows the 
XIS-0 image of Cen A before and after the pointing-error correction. The image is sharpened 
by the pointing-error correction. 

4. Point Spread Function Modeling 

Four XRTs have different complex position-dependent PSFs and they were calibrated on 
ground (Itoh et al. 2004a; Itoh et al. 2004b; Misaki et al. 2004). The PSFs observed in orbit are 
different from those measured on the ground because of absence of gravity in orbit. Hence, the 
PSF models made on the ground calibration have to be modified for the present application. 

We constructed two PSF models: one based on the ray-tracing simulator, xissim, which 
we call the xissim PSF and the other on the observed Cen A image, the observed PSF. The 
xissim PSF implements the position-dependence in the entire image area. However xissim is 
based on the ground calibration and does not agree with the in-orbit point source images. 
On the other hand, the observed PSF is available only near the optical axis of each XIS. We 
combined the two PSF models when possible: the observed PSF for the central region within 
6' from the optical axis and the xissim PSF for the outer region. 

6 



-20 















i 


— i h 



10 20 

Time [hr] 



Fig. 1. Time variation of the position of Cen A in the XIS-0 image during ob- 
servation started at 2005-08-19 05:55:00 (UTC). The CCD pixel width is 1".04. 
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Fig. 2. XIS-0 images of Cen A in 0.5-10 keV band before the 
XRT alignment-error correction (left) and after the correction (right). 

4.1. XissimPSF 

The XRT ray-tracing library, xissim, is included in Suzaku FTOOLS in the HEADAS 
software package (Ishisaki et al. 2007). We simulated PSF images for 16 x 16 source locations, 

(DETX, DETY) = (32 + 64 x z, 32 + 64 x j) (0 < i < 16, < j < 16) (13) 

where each PSF represents a l'.l x 1M segment on the 17' .7 x 17' .7 XIS image area. We 
assumed that the four quadrants of each XRT mirror are identical and each quadrant has a 
mirror-symmetry relative to the median angle of the quadrant. This means that the position 
dependence of PSF has four-fold axis-symmetry around the optical axis of the XRT and each 
quadrant has internal mirror symmetry. The position dependence of PSF has been simulated in 
the 1/8 triangular region of the square XIS FOV. The left panel of figure 3 shows the simulated 
position dependence of PSF. 

The simulated PSF images have then been fitted with an analytic model function. Figure 
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Fig. 3. (Left): PSF samples generated by XRT/XIS simulator, xissim. Open triangles are 
positions at which PSFs are simulated. Open circles show the grid points where PSF can 
be taken from the symmetric point in the triangle. (Right): PSFs generated with our 
Xissim PSF model. Gray contour colors are spaced logarithmically in the both figures. 

4 illustrates parameters used in the PSF model: the X-ray source position (0,0), the detected 
photon position (p, < 0), and coordinates for the mirror quadrant boundaries by (a, (3) The 
function of the xissim PSF, PSF^^O^] PjVO? consists of a central peak profile represented 
by three exponential distributions, p(6,(f)] p,V0- The sensitivity gaps at the boundaries of the 
mirror quadrant are represented by q(9,(j)] p,ip). They take the following formulae: 

PSF^ssimiO, (j)] p, ip) = p(0, (j)] p, ip)q(0, (j)] a, (5) 



p(0,0;p,^) 




g(0,0;a,/3) = 



1 + exp 



1 + exp 



of cos 2 ip + 6? sin 2 
k 2 {\a\-k 1 )-\(5\ 



(14) 
(15) 



^ fc 3 (|a| - fci) 
'Wl-fci) 




(M<*i, I 

(H>l/?l) 
(M<lfl) 



(16) 



Each simulated PSF at 



\h(\P\-h) + h A 

is fitted to obtain the 17 parameters, {a^, fo^, Cj}i=i,2,3 an d 
{fci, ^i}i=i,2,3,4- The right panel of figure 3 shows the fitted PSF. We calculate the PSF for 
arbitrary position by interpolating from the simulated PSFs at the grid points. 

4-2. Validation of xissim PSF 

We have validated the xissim PSF model by comparing actual observed image of Cen 
A in 2-10 keV band. The central X-ray source is known to be point-like in the energy band 
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Fig. 4. Coordinate definition used in PSF modeling. 



(Evans et al. 2004; Markowitz et al. 2007). 

Figure 5 shows the Cen A images taken with the four XISs after the XRT alignment-error 
correction, and the model PSF. The observed Cen A image is similar to the PSF model, but 
there are apparent differences. Figure 6 shows the radial profiles of the observed Cen A images 
and the xissim PSF model averaged over the 360° azimuthal angle and the azimuthal profiles of 
an annulus between radii r = 40 pixel and r = 60 pixel. The azimuthally averaged radial profile 
of the Cen A images agree well with the model PSF except for the central core of radius r < 5". 
On the other hand, the azimuthal profiles of the Cen A images are significantly different from 
the model PSF, reflecting the complicated asymmetric profiles of the real XRT-XIS system. 

Among the Cen A images observed by the four XISs, that of the XIS-1 came closest to 
the model PSF. Hence, we have decided to use XIS-1 to test the image deconvolution method 
in section 5. 

4.3. Observed PSF 

We made another PSF model out of the observed Cen A image, and call it the observed 
PSF. We plot radial profiles of the Cen A image for several azimuth angles in the left panel of 
figure 7. These profiles are fitted with a model functions consisting of two exponentials. The 
PSF model is formulated on the best-fit radial profiles. The obtained observed PSF model is 
shown in the right panel of figure 7. 

The archived data of Cen A were taken in one pointing and hence position dependence 
of PSF could not be extracted. Because of this, we use the observed PSF only for the central 
region of the XIS FOV in the deconvolution analyses described below. 
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Fig. 6. (Left): Radial profiles of xissim PSF and images of Cen A in 2-10 keV ob- 
served by each XIS detector. All profiles are normalized to 0.01 at a radius of 60 pixel. 
(Right): Azimuthal profiles extracted from annulus between radii r = 40" and r = 60". 
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Fig. 7. (Left): Radial profiles of XIS-1 Cen A image in 2-10 keV band with 
statistical errors and their best-fit model functions (solid lines) for various az- 
imuthal angles. (Right): Observed PSF model for central region of XIS-1 image. 



5. Application to XIS images 

We have applied the present deconvolution method to simulated XIS images as well as 
to observed images in the Suzaku archival data. We have chosen Cen A, PSR B1509-58 and 
RCW 89 as targets because they include both a bright point source and extended emissions 
and also because there are Chandra observations of the same regions. The bright point sources 
serve to correct for the XRT pointing error as has been described in section 3.1. The images of 
Chandra ACIS are used to evaluate the fidelity of the deconvolved images. 

5.1. Deconvolution of Images 

The deconvolution program is developed on a standard Linux machine in the C program 
language using BLAS 1 / LAPACK 2 linear algebra program libraries (Anderson et al. 1999) 
customized by ATLAS project (Whaley et al. 2001). 

The deconvolution of XIS images proceeds in the following steps. 

1. Extract a region of interest from the entire XIS image (1024x1024 pixels covering area of 
17'. 8 x 17'. 8) and make a 64x64 tiled image where each tile is combination of 6 x 6 XIS 
CCD pixels (each CCD pixel covers 6". 3 x 6". 3). For RCW 89, the tile size has been set 
larger to cover 8x8 CCD pixels because number of photon per tile was lower. 

2. Make a response function based on the two PSF models. The PSF is normalized to 1.0 so 

1 http://www.netlib.org/blas/ 

2 http: / / www.netlib.org/lapack/ 
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that the number of photons is preserved. We employed the observed PSF model for the 
region within 6' from the XRT optical axis. Exception is RCW 89. The xissim PSF model 
was used for RCW 89. 

3. Calculate the inverse of the response matrix. 

4. Multiply the inverse response matrix with the tiled image 

5. Smooth the response-inverted image adaptively. The adaptive smoothing takes a signal- 
to-noise ratio, SNR opt as a parameter. We set SNR opt =4.0 as described in section 2.3 for 
all the examples described below. 

5.2. Test with Simulated Images 

We have tested the fidelity of the present deconvolution method using simulated XIS-1 
images. The simulated image consists of three point sources and an extended source as shown 
in figure 8. The surface brightness of the extended emission is represented by a simple j3 model: 



where the critical radius is set to r c = 30" and (5 to 0.5. Relative fluxes of three point sources and 
extended emission are set at 4 : 2 : 1 : 40. To study how fidelity depends on photon statistics, we 
have deconvolved three images containing total signal photons of 20,000, 100,000, and 500,000 
and background photons expected for nominal blank-sky in a 50-ks exposure (Mitsuda et al. 
2007). These images are then compared with images deconvolved by a Richardson-Lucy method 
with 100 iterations in figure 8. 

One can see that the Richardson-Lucy method adds photons to high points in Poisson 
fluctuation noise and depletes photons from low points in the fluctuation. Our method, on 
the other hand, smooths out the fluctuation below the predetermined signal-to-noise ratio 
while reproducing the three point sources and the extended source well. Such noise filter can 
be added to the Richardson-Lucy method to suppress artifacts but the filtering strategy is 
strongly coupled with the number of iteration. We note that the raw image of 100,000 photons 
at the middle row in Figure 8 has about ^ 50 signal photons and ~ 3 noise photons per tile in 
average within the critical radius r c = 30". 

5.3. Cen A 

The data of Cen A are divided into two energy bands, > 3 keV and < 3 keV, and the 
two images are deconvolved separately. In the hard band, the flux of Cen A is known to be 
dominated by the central point source (Evans et al. 2004; Markowitz et al. 2007). In the soft 
band, an extended jet profile with a scale of 30^-lSO" has been observed by Chandra (Kraft et 
al. 2002; Kataoka et al. 2006). 
5.3.1. Hard band: 3-10 keV 

Figure 9 shows the observed image of XIS-1 in 3-10 keV band, its deconvolved image, 
and the Chandra ACIS image rebinned to the tile size. We have also convolved the Chandra 




(17) 



12 




Fig. 8. Deconvolution of simulated XIS-1 images. (Left) Model image consists of three point 
sources and an extended emission. Relative fluxes of the three point sources (A, B, C) and ex- 
tended emission are 4:2:1: 40. (Center left) Simulated XIS-1 images with photon counts of 
20,000, 100,000, 500,000 from top to bottom. Background photons with Poisson distribution ex- 
pected for a nominal blank-sky in a 50-ks exposure are included. (Center right) Deconvolved 
XIS-1 images with 20,000, 100,000, and 500,000 photons from top to bottom by the present 
method. (Right) Deconvolved XIS-1 images by a Richardson-Lucy method after 100 iterations. 
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image with the XIS-1 PSF and shown in the top right panel. We note that the Chandra ACIS 
image is distorted by the pile-up effect at the Cen-A nucleus and the one shown in figure 9 has 
been corrected for the distortion by using the Suzaku XIS image. 

We have studied the deconvolved Suzaku image by slicing it along the 3 blue dash-dot 
lines drawn at (f) = — 33°, 0°, and 33° from the Right Ascension axis crossing the Cen A core in 
figure 9. The surface brightness profiles are compared with those of the XIS-1 raw image and 
the Chandra ACIS image in figure 10. The surface brightness was calculated by assuming a 
power-law spectrum with photon index Y = 1.8 (Markowitz et al. 2007). The vertical scale of 
the XIS-1 raw image was normalized at the peak to the deconvolved image. 

In the deconvolved image, the large wings extending over ~ 120" in the raw images 
are drastically reduced. The width of the central peak in the deconvolved image is wider 
by 12" — 24" than the ACIS image. This is interpreted that the spatial resolution has been 
improved to 12" — 24" by deconvolution but not better. The spatial resolution depends on the 
azimuth angle: It is the worst at (j) = 33°, which corresponds to the direction along which the 
PSF is extended. We find two artifact peaks in the region 50" — 100" away from the peak in 
all three profiles at about 1/50 of the peak surface brightness. This shows limitation of the 
present deconvolution method. With the present modeling of the XIS PSF, our deconvolution 
reconstructs the image down to about 1:50. A third exponential component with a wider wing 
may be required to improve the dynamic range. 
5.3.2. Soft band: 0.5-3 keV 

Figure 11 shows the observed image of XIS-1 in 0.5-3 keV band, the deconvolved image, 
and the rebinned Chandra ACIS image. The convolved ACIS image with the XIS-1 PSF is also 
shown in the top right panel. The pile-up effect at the nucleus of Cen A in the ACIS image is 
corrected from the flux measured by the XIS. Our deconvolution reconstructs the bright blob 
a B" but not the blob a C" in the jet. 

Figure 12 shows the cross-section profiles of the XIS-1 raw, deconvolved, and ACIS 
images sliced along the direction to the north-east jet (a dash-dot blue line in figure 11). Here, 
the surface brightness is calculated assuming a power-law spectrum with photon index T = 1.3 
(Markowitz et al. 2007). The profile around peak is well restored in the deconvolved image 
with a spatial resolution of 20". The restored image is consistent with that by Chandra if we 
take the pile-up effect in the Chandra image into account. 

5.1 PSR B1509-58 and RCW 89 

The observed region of PSR B1509-58 includes a bright point source PSR B1509-58 and 
an extended source RCW 89 (Yatsu et al. 2005; DeLaney et al. 2006). Figure 13 shows the 
overall image by XIS-1. The hard band image exhibits only one point source, PSR B1509-58, 
and no extended emission. The image in the soft band consists of complex extended emissions 
as shown in the next subsection. 
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Fig. 9. XIS-1 raw image (top left), convolved Chandra ACIS image with XIS-1 PSF (top right), 
XIS-1 deconvolved image (bottom left), and Chandra ACIS image (bottom right) of Cen A in 3-10 
keV. All images are binned with a same unit tile size of ~ 6" x 6". Contour colors are spaced 
logarithmically. Notice that the Chandra ACIS image suffers from the pile-up effect at the nu- 
cleus of Cen A (the brightest point). Blue dash-dot lines indicate the sliced directions in figure 10. 




arc-second arc-second arc-second 



Fig. 10. Cross-section profiles of the XIS-1 raw image (red dot) deconvolved image (black solid) 
and Chandra ACIS image (blue dash) of Cen A in 3-10 keV band, sliced along the blue dash- 
dot lines in figure 9. Error bars represent 1-a photon-statistics uncertainties. The vertical scale 
of the raw image is normalized by the peak value to the deconvolved image. Notice that the 
pile-up effect at the peak of the Chandra ACIS image (one pixel at the center) is corrected. 
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Fig. 11. XIS-1 raw image (top left), convolved Chandra ACIS image with XIS-1 PSF (top right), 
XIS-1 deconvolved image (bottom left), and Chandra ACIS image (bottom right) of Cen A in 0.5-3 
keV band. All images are binned with a same unit tile size of ~ 6" x 6". In images of XIS-1, 
the blank areas at the top left and the bottom right are not covered by clocked CCD window in 
this observation. Note that the peak of the Chandra ACIS image (one pixel at the center) suf- 
fering from a pile- up effect has been corrected using the flux measured by the Suzaku XIS. The 
blue dash-dot lines on the XIS-1 deconvolved image and the Chandra image indicate the sliced 
direction in Figure 12. The labels "A", "B", "C" in the Chandra image identify the peaks. 

5.4.1. PSR B1509-58: 0.5-3 keV 

The raw image of the region including PSR B1509-58 and RCW 89 is shown in figure 
13. The image includes interesting extended emissions and a point source, and can be viewed 
as a typical target for our image deconvolution program. 

Figure 14 shows the raw XIS-1 image and the deconvolved image in the 0.5-3 keV band, 
and the Chandra ACIS image with that convolved with the XIS-1 PSF. These cross-section 
profiles along the Right Ascension axis are also shown. The surface brightness was calculated 
assuming a power-law spectrum with photon index T = 1.8. The Chandra image suffers from 
the pile-up effect at the position of PSR B1509-58. The sliced surface brightness profile of the 
deconvolved image is consistent with the Chandra image if the pixel saturation is corrected 
using the XIS-1 image. 
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Fig. 12. Cross-section profiles of XIS-1 raw image (red dot) and deconvolved image (black 
solid) and Chandra ACIS image (blue dash) of Cen A in 0.5-3 keV sliced along the direc- 
tion of the jet (blue dash-dot lines in figure 11). Error bars represent 1-a photon-statis- 
tics uncertainties. The vertical scale of XIS-1 raw data is normalized at the peak value to 
the deconvolved image. The peak of the Cen- A nucleus in the Chandra profile has been cor- 
rected for the pile- up effect. The labels "A", "B", "C" indicate the peaks in figure 11. 




Fig. 13. XIS-1 overall raw image of PSR B1509-58 and RCW 89 in 0.5-3 keV band. 
Solid boxes represent the regions of which images are applied to the deconvolution analysis. 
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Fig. 14. PSR B1509-58 region: XIS-1 raw image (top left), convolved Chandra ACIS image with 
XIS-1 PSF (top center), XIS-1 0.5-3 keV deconvolved image (bottom left), Chandra ACIS 0.5-3 keV 
image (bottom center), and cross-section profiles sliced along Right Ascension (blue dash-dot lines in 
the images) (bottom right). Both XIS-1 and Chandra ACIS images are binned with a same tile size 
of ~ 6" x 6". A pixel of the central core of PSR B1509-58 in the Chandra image is corrected for 
the pile-up effect. Error bars in cross-section profiles represent 1-a photon-statistics uncertainties. 

The jet-like structure extending from the pulsar is restored in the deconvolved image. 
However other extended features around the pulsar are missed or incorrectly restored. Their 
surface brightness is lower than 1/50 of the peak value, or below the dynamic range of this 
deconvolution method. 
5.4.2. RCW 89: 0.5-3 keV 

Figure 15 shows the raw XIS-1 image and the deconvolved image of RCW 89 in the 
0.5-3 keV band, and the Chandra ACIS image with that convolved with the XIS-1 PSF. These 
cross-section profiles along the Right Ascension axis are also shown. For this source, 8x8 
raw pixels were combined to one tile to secure high numbers of photons per tile. The xissim 
PSF model was used in the response matrix inversion because the image area extended larger 
than 6' from the XRT optical axis. The surface brightness is calculated assuming a power-law 
spectrum with photon index Y = 1.8. The extended bright region of scale greater than ~ 20" 
in the upper part of the image has been reproduced well in the deconvolved XIS-1 image but 
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Fig. 15. RCW 89 region: XIS-1 raw image (top left), convolved Chandra ACIS image with XIS-1 
PSF (top center), XIS-1 0.5-3 keV deconvolved image (bottom left), Chandra ACIS 0.5-3 keV im- 
age (bottom center), and cross-section profiles sliced along Right Ascension (blue dash-dot lines in 
the images) (bottom right). Both XIS-1 and Chandra images are binned with a same tile size of 
8" x 8". Error bars in peak cross-section profiles represent 1-cr photon-statistics uncertainties. 



positions of narrower high points do not agree with those in the Chandra image. We interpret 
this due both to the inaccuracy of the xissim PSF model and to the limited angular resolution 
of ~ 20". Within these limitations, the deconvolution method has reproduced the prominent 
shell-like structure in the upper part. 

6. Conclusion and Future Prospect 

We have developed an image deconvolution method for the Suzaku XIS based on response 
matrix inversion and adaptive smoothing. The method has been tested with two XIS-1 images 
both containing extended sources and one prominent point source: Cen A and PSR B1509- 
58/RCW 89. By comparing the deconvolved images with the corresponding Chandra ACIS 
images, we conclude that spatial resolution has been restored to ~ 20" to a brightness level 
around 1/50 of the brightest tile in the image. 

Recent X-ray instruments including Suzaku are finding complex morphology of thermal, 
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non-thermal and K-shell line emissions in many extended sources including young supernova 
remnants (SNRs) (e.g. Ueno et al. 2007), pulsar wind nebulae (PWNs) (e.g. Seward et al. 
2006), the Galactic Center region (e.g. Koyama et al. 2007b; Koyama et al. 2007c; Koyama 
et al. 2007d), and galaxy clusters (e.g. Sanders et al. 2005). To understand such objects, the 
spatial resolution of Suzaku XIS has to be improved substantially. The present work provides 
a procedure for such improvement. 

We plan to improve PSF modeling and incorporate the XRT alignment scheme developed 
by the Suzaku XIS team. 
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